#include "stdafx.h"

#include "Solver.h"
#include <math.h>
#include "matrix.hpp"
#include "evector.hpp"
Solver::Solver()
{
    //ctor
}

Solver::~Solver()
{
    //dtor
}


EVector<double> Solver::CGMethod(FullMatrix<double> B, EVector<double> F, double tol, int * nItr)
{
	int degree = F.size();
    EVector<double> _x(degree);
    EVector<double> _r(degree);
    EVector<double> _p(degree);
    EVector<double> _w(degree);
    EVector<double> temp(degree);
    _r = F;
    double rho0, rho1;
    double alpha;

	rho1 = _r.GetNorm2();
    rho0 = 0;
    int count = 1;
    while(sqrt(rho1) <= tol*F.GetNorm2())
    {
        if(count == 1)
            _p = _r;
        else
            _p = (_r + _p * (rho1/rho0));

        _w = (B * _p);
        alpha = rho1/_p.ScalarProd(_w);

        _x=(_x + _p * alpha);
        _r= (_r - _w * alpha);

        rho0 = rho1;
        rho1 = _r.GetNorm2();
        count++;
    }
    *nItr = count;
    return _x;
}
